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EXECUTIVE  SUMMARY 


In  this  research  we  have  successfully  carried  out  an  implementation  of  the  differen¬ 
tial  inversion  algorithm  developed  by  Dr.  Jean  I.  F.  King  of  the  Air  Force  Geophysics 
Laboratory.  This  algorithm  is  an  exciting  advance  in  the  methodology  of  atmospheric 
temperature  sensing  because  it  allows  the  computation  of  atmospheric  temperature  pro¬ 
files  from  upwelling  radiances  and  their  derivatives  in  a  computationally  effective  fashion. 

The  primary  problem  in  any  implementation  of  the  differential  inversion  is  the  control 
of  numerical  instability.  A  portion  of  this  instability  is  intrinsic  to  the  differential  inversion 
algorithm  because  it  is  based  upon  Laplace  transform  techniques;  the  rest  is  due  to  use 
of  finite  difference  schemes  to  evaluate  numerical  derivatives  (which  are  sensitive  to  the 
presence  of  experimental  and  other  errors) .  We  have  learned  how  to  control  these  numerical 
instabilities  as  a  result  of  studies  carried  out  on  synthetic  radiance  data  which  characterized 
the  properties  of  the  differential  inversion  algorithm  on  a  controlled  data  set. 

We  then  applied  the  differential  algorithm  to  upwelling  radiance  data  from  SCRIBE 
and  TOVS.  Owing  to  the  difficulties  in  calibration  and  the  noise  level  in  the  SCRIBE  data, 
inversion  errors  on  the  order  of  20  K  were  obtained.  However,  inversion  of  TOVS  data 
involved  errors  on  the  order  of  only  2  K. 

We  have  demonstrated  that  the  differential  inversion  algorithm  may  be  successfully 
employed  to  invert  radiance  data  to  obtain  atmospheric  temperatures.  Our  work  suggests 
several  profitable  avenues  for  research  to  improve  the  accuracy  and  utility  of  the  differen¬ 
tial  inversion  algorithm,  both  through  improvements  in  the  algorithm  itself  and  throught 
developments  in  temperature  sounder  instrumentation. 


1  INTRODUCTION 


Satellite  observations  provide  the  only  practical  means  of  obtaining  temperature  pro¬ 
files  in  the  atmosphere  over  large  geographic  areas  (especially  over  oceans)  and  over 
long  periods  of  time.  This  implies  that  the  development  of  new  and  computationally 
effective  algorithms  for  measuring  atmospheric  temperature  profiles  from  satellite  data 
is  a  significant  and  important  task. 

The  measurable  quantity  most  accessible  to  satellite  observations  is  the  infrared 
upwelling  radiance.  Most  typically  this  is  measured  at  wavelength  ~  15  ^m,  near 
the  atmospheric  COj  absorption  band.  Different  spectral  regions  in  this  wavelength 
region  become  optically  thick  at  different  heights  in  the  atmosphere,  and  so  the  up- 
welling  infrared  radiance  in  those  wavelength  intervals  contains,  at  least  in  principle, 
information  about  the  Planck  function  at  different  heights  in  the  atmosphere.  On  this 
basis,  one  should  be  able  to  determine  a  temperature  profile  of  the  atmosphere  using 
measurements  of  the  upwelling  infrared  radiance. 

In  practice,  we  must  solve  the  radiative  transfer  equation  for  the  Planck  function 
over  this  spectral  interval,  a  Fredholm  integral  equation  of  the  first  kind.  Since  we  are 
dealing  with  an  integral  equation,  we  must  cope  with  problems  of  non-uniqueness  of 
the  solution  obtained.  This  is  a  problem  of  practical  importance  because  the  range  of 
plausible  solutions  of  the  transfer  equation  (within  limits  imposed  by  the  accuracy  of 
the  experimental  data)  cover  an  unacceptably  large  range  of  temperatures. 

A  number  of  workers  have  approached  this  problem  with  techniques  designed 
to  obtain  the  “true”  temperature  profile  (or  at  least  a  pragmatically  useful  profile) 
by  a  variety  of  methods,  most  of  which  require  some  a  priori  assumptions  such  as 
the  assumption  of  an  initial  temperature  profile.  Application  of  relaxation  techniques 
(Chahine  1970,  1972;  Smith  1970;  Yeh,  Vonder  Haar,  and  Liou  1985)  requires  the 
selection  of  an  initial  temperature  profile  to  commence  an  iterative  process.  Application 
of  the  method  of  Backus  and  Gilbert  (Conrath  1972)  requires  the  assumption  of  a 
reference  temperature  profile  in  order  to  evaluate  radiative  transfer  kernels. 

Other  relevant  work  has  been  undertaken  (Twomey  1970)  to  delineate  the  limits 
in  information  content  of  solutions  of  the  radiative  transfer  equation.  In  this  vein, 
Fourier  analysis  of  the  radiative  transfer  equation  (Gautier  and  Revah  1975)  allows 
treatment  of  the  temperature  profile  as  a  band-limited  spectrum. 

In  this  paper  we  consider  a  new  approach  to  this  problem  based  on  Laplace 
transform  theory  (King  1983,  1985).  This  technique  uses  the  upwelling  radiance  mea¬ 
surements  and  their  logarithmic  pressure  derivatives  to  construct  an  atmospheric  tem¬ 
perature  profile.  This  method  has  the  particular  theoretical  attraction  of  assuming  no 
a  priori  information  about  the  temperature  profile  beyond  the  fact  that  it  is  reason- 


ably  smooth  (i.e.  analytic).  Our  numerical  studies  and  simulations  using  this  algorithm 
have  investigated  the  viability  of  this  approach  and  indicated  some  of  its  sensitivities  to 
truncation  error,  roundoff  error,  experimental  uncertainty,  and  other  noise  processes. 
This  differential  inversion  technique  has  been  applied  to  a  set  of  Fourier  transform 
spectrometer  data  in  the  15  (j, m  spectral  region  and  to  a  sample  of  data  from  the 
NOAA  TIROS  Operational  Vertical  Sounder  (TO VS).  On  the  basis  of  these  studies, 
we  shall  evaluate  this  technique  and  discuss  instrument  designs  optimally  suited  to  its 
application. 


2  THEORY  OF  DIFFERENTIAL  INVERSION 


2.1  General  Theory 

In  this  section  we  develop  the  theory  of  differential  inversion  techniques  as  applied 
to  the  detemination  of  atmospheric  temperature  profiles  on  the  basis  of  the  infrared 
upwelling  radiance.  This  infrared  upwelling  radiance  might  be  measured  by  a  satellite 
or  a  balloon-borne  spectrometer  at  high  altitudes.  For  sake  of  completeness,  we  shall 
develop  the  theory  of  differential  inversion  in  some  detail.  Our  treatment  follows  the 
development  of  this  subject  by  King  (1985). 

Radiative  transfer  theory  relates  the  upwelling  thermal  radiance  of  an  atmo¬ 
sphere  to  an  integral  transform  over  pressure  of  the  Planck  intensity  B(p)  and  a  weight 
function  ^(p/p).  The  weight  function  is  determined  by  the  detailed  line  character  of 
the  sampled  interval,  effectively  sampling  that  height  of  the  atmosphere  at  which  the 
optical  depth  in  that  spectral  interval  is  near  unity.  Therefore,  the  radiance  of  the 
wavelength  channel  whose  weight  function  peaks  at  p  =  p  is  given  by 

£(?)  =  /  S(p)W(p/p)dp/p.  (1) 

Jo 


A  change  of  independent  variable  into  logarithmic  pressure  units,  with  modi¬ 
fied  functions  B(p)  =  5(-lnp),  R(p)  =  J?(-lnp),  and  W(p)  —  iy(lnp)  converts  the 
expression  for  the  upwelling  radiance  (1),  into  a  convolution  integral, 

R[s)  =  r°  B[z)W{$  -  z)  dz  =  B  *W,  (2) 

J  —  OO 

where  z  —  —  lnp,  and  f  =  —  lnp  are  height-like  variables.  This  convolution  integral 
form  motivates  the  application  of  Laplace  transform  theory  to  solve  for  the  Planck 
function  and  equivalently  for  the  temperature.  Application  of  the  bilateral  Laplace 
transform  using  the  Faltung  theorem  results  in 

r(s)  =  b(s)w(-  s),  (3) 


with  r,  6,  and  w,  the  Laplace  transforms  of  R,  B ,  and  Wy  respectively,  and  s  the 
Laplace  transform  variable.  It  will  be  seen  that  the  bilateral  transform  of  the  weight 


function 


tu(-s)  =  f  e  ,zW(z)  dz  —  f  p  ’W(p)  dp/p  ==  0(1  —  s), 

J  -  oo  J  0 


may  be  viewed  as  a  generalized  T  function.  This  is  significant  for  practical  applications 
of  the  theory  in  which  it  may  be  desirable  to  approximate  atmospheric  weight  functions 
in  some  convenient  fashion. 


We  shall  use  expansion  techniques  to  solve  for  the  Planck  function  as  a  function 
of  height.  Proceed  by  expanding  the  reciprocal  of  the  fl  function  in  a  MacLaurin  series, 

n(lbj  =  S^>‘ 

in  which  the  coefficients  A*  are  given  by 

A‘  -  — ki — '  <6> 

The  expansion  coefficients  are  thus  seen  to  involve  the  A:th  derivative  with  respect  to 
the  transform  variable  s  of  the  bilateral  weight  transform,  evaluated  at  s  =  0.  We  may 
now  express  the  Laplace  transform  of  the  Planck  function  in  terms  of  a  series  involving 
the  Afc’s  and  the  transform  of  the  radiance, 

b{s)  =  f^\kskr{s).  (7) 


Taking  the  inverse  Laplace  transform  of  (7)  then  yields  the  differential  inversion  formula 
we  seek, 

B(z)  =  f;  A *£><*>  12(<r)  -  J2(f)  +  A +  A 3JT(c)  +  ...  (8) 

t= o 

This  formula  has  been  derived  formally  for  Laplace  transform  inversion  by  Widder 
(1971). 


2.2  Differential  Inversion  Coefficients  and  Generalized  Weight 
Functions 

Calculation  of  the  Planck  function  profile  depends  on  the  inversion  coefficients  which 
are  determined  by  the  particular  form  of  the  weight  function.  Determination  of  de¬ 
tailed  atmospheric  weight  functions  over  a  range  of  frequencies  is  a  computationally 
demanding  task.  It  is  frequently  desirable  to  work  with  an  approximate  form  which  has 
most  of  the  relevant  features  of  actual  atmospheric  weight  functions  and  the  advantage 
of  computational  tractability. 

We  shall  now  develop  closed-form  expressions  and  numerical  values  using  such 
an  approximation  to  the  weight  functions.  Define  the  generalized  exponential  weight 
function, 

Vm{P)  =  T,  +  exP{-mP'^),  (0) 

where  P  —  p/p,  with  p  the  pressure  at  which  the  weight  function  achieves  its  maximum. 
The  generalized  exponential  weight  function  satisfies  the  normalization  requirements, 


The  parameter  m  adjusts  the  effective  width  of  the  generalized  exponential  weight 
function,  smaller  values  of  m,  “sharpen”  the  weight  function,  t.e.  decrease  its  effective 
width.  It  should  be  noted  that  the  particular  choices  of  m  =  1  and  m  =  1/2  correspond 
to  the  Goody-statistical  and  the  Elsasser  band  models,  respectively. 


The  bilateral  Laplace  transform  of  the  weight  function  follows  (see  eq  (4)  )  as 
f°°  e—Wm{z)  dz  =  r  P- Wm(p)  dp/p  =  rm(l  -  s),  (1 

J  —  oo  JO 

where  we  have  defined  the  mth  order  T  function  as 


rm(s  + 1) 


_  r[m(s  -I-  1)] 


The  evaluation  of  the  coefficients  A*,  as  can  be  seen  from  Equations  (4),  (6),  and 
(12),  involves  the  derivatives  of  the  T  function.  The  first  derivative  is  related  to  the 
digamma  or  xj}  function 

*W  =  r'(s)/r(s),  (is) 

and  its  higher  derivatives  the  polygamma  functions 

1 

^(s)  =  -~~-\nT{s).  (14) 

Closed  form  expressions  for  the  inversion  coefficients  A k  (k  =  0, 1, ...  ,6)  as  a  function 
of  the  weight  parameter  m  are  given  in  King  (1985).  For  integer  and  half-integer 
argument,  the  polygamma  functions  are  expressible  in  terms  of  the  Riemann  zeta 
functions. 

Calculation  of  inversion  coefficients,  given  for  some  special  cases  in  Tables  1  and  2, 
is  straightforward  but  tedious.  Note  that  as  the  weight  function  becomes  increasingly 
concentrated,  corresponding  to  smaller  m,  that  all  A*’s  tend  to  zero  for  k  >  1,  while 
Ai(m  — ►  0)  — >  —1.  This  implies  that  even  for  a  delta-type  weight  function,  there  is  a 
shift  between  the  Planck  intensity  and  the  radiance  profile  extrema,  t.e.  pmin  pm,„, 
where  Bmin  =  B(pmm),  and  £m,„  =  ^(pn,„).  Furthermore,  the  extent  to  which  the 
odd  coefficients  differ  from  zero  is  a  measure  of  the  asymmetry  of  the  weight  function. 
Therefore, 

A,^0  (t  =  l,3,5,...).>W'(-2)#iy(z). 


2.3  Expected  Limitations  of  Differential  Inversion  Techniques 

This  theoretical  discussion  of  the  differentia!  inversion  technique  allows  us  to  antici¬ 
pate  some  of  the  difficulties  that  might  be  encountered  in  implementing  this  method 


Table  1:  Inversion  coefficients  for  generalized  exponential  weight  function  with  m  =  1 
taken  from  (King  1985). 


+  1.0 

-0.5772150649 

-0.6558780715 

+0.0420026350 

+0.1065386114 

+0.0421977346 

-0.0096219715 

-0.0072189432 

-0.0011651076 

+0.0002152417 

+0.0001280503 

+0.0000201349 

-0.0000012505 


Table  2:  Differential  inversion  coefficients  for  representative  general  exponential  weight 

functions. 

Delta 

Strongly 

Elsasser 

Statistical 

Strong 

function  peaked 

clustering 

m  =  0 

m  =  0.1 

m  =  0.5 

m  =  1 

m  =  4 

Ao 

+  1.0 

+  1.0 

+  1.0 

+  1.0 

+  1.0 

Ai 

q 

1 

-0.8121169848 

-0.6351814227 

-0.5772156649 

-0.5207067709 

a2 

0.0 

-0.1773994973 

-0.4151225552 

-0.6558780715 

-2.1350158750 

A3 

0.0 

-0.9119348444 

-0.0014993282 

+0.0420026350 

0.3050207816 

A4 

0.0 

+0.0004253445 

+0.0416237314 

+0.1665386114 

2.2390138114 

A5 

0.0 

+0.0001167500 

+0.0104035581 

+0.0421977346 

+  0.4628141296 

for  obtaining  atmospheric  temperature  profiles.  First  of  all,  it  is  apparent  that  we  will 
obtain  some  level  of  truncation  error  owing  to  the  necessity  of  terminating  the  series 
in  Equation  (8)  after  some  finite  number  of  terms.  In  practice,  the  number  of  terms 
at  which  we  will  have  to  terminate  this  series  will  depend  on  the  reliability  of  our 
measured  radiances  and  their  derivatives,  and  as  is  well  known,  numerical  differenti¬ 
ation  of  experimental  data  corrupted  by  noise  is  at  best  a  chancy  business.  It  would 
be  very  surprising  if  more  than  the  first  few  derivatives  of  radiance  data  from  actual 
measurements  were  reliable.  Therefore,  we  will  investigate  the  extent  to  which  this  is 
of  importance  in  the  differential  inversion  algorithm. 

Practical  difficulties  also  arise  in  the  determination  of  the  weight  functions.  The 
weight  functions  may  either  be  determined  by  regression  techniques  or  determined  by 
differentiation  of  atmospheric  transmittance  profiles.  Therefore,  there  is  some  uncer¬ 
tainty  associated  with  the  weight  function,  and  with  the  location  of  its  maximum,  p. 
This  is  especially  true  when  the  weight  function  is  obtained  by  differentiation  of  the 
transmittance,  where  problems  in  accurately  obtaining  numerical  derivatives  may  be 
important.  The  uncertainty  in  the  weight  functions  is  also  reflected  in  the  evaluation 
of  the  A  coefficients  which  depend  on  moments  of  the  weight  function. 

In  addition  to  these  practical  limitations,  a  more  fundamental  difficulty  must  be 
considered.  Evaluation  of  inverse  Laplace  transforms,  as  is  required  in  differential  in¬ 
version,  is  known  to  be  a  numerically  unstable  process  (Widder  1971),  compared  to 
the  situation  in  Fourier  inversion,  where  demonstrably  stable  and  accurate  algorithms 
exist.  The  reason  for  this  numerical  instability  arises  ultimately  from  the  fundamental 
different  character  of  the  exponential  function  for  real  and  imaginary  arguments.  For 
real  arguments,  z,  e~z  is  a  smooth  and  monotonic  function  from  x  —  0  to  x  =  oo; 
and  for  imaginary  arguments  the  exponential  becomes  oscillatory,  maintaining  its  av¬ 
erage  amplitude  for  all  z.  It  is  apparent  then  if  we  want  the  Laplace  transform, 
/0°°  e~tzf(x)  dx,  to  show  some  response  to  values  of  z  around  z  =  100,  say,  then  s  must 
be  small  enough  that  e~,z  is  appreciably  different  from  zero.  Say  that  we  need  $  «  0.01; 
now  the  derivative  of  e~,z  is  —  se~,z,  and  since  s  is  small,  then  so  is  its  derivative  and 
so  the  function  e~,z  changes  only  slowly  around  z  —  100.  This  implies  that  we  will 
have  considerable  difficulty  in  resolving  features  in  the  function  /(z)  around  z  =  100, 
i.e.  distinguishing  /( 100)  from  say  /(99)  or  /(101).  Also  rapid  fluctuations  in  /(z)  in 
the  neighborhood  of  z  =  100  will  give  only  a  small  change  in  J£°  e~,z  f(x)  dx.  These 
considerations  show  that  kernels  with  smooth  behavior  can  yield  unstable  solutions, 
i.e.  in  terms  of  the  physical  variables  of  our  particular  problem,  small  changes  in  the 
upwelling  radiance  can  yield  large  changes  in  the  inferred  temperature  profile  one  ob¬ 
tains  as  a  result  of  differential  inversion  computations.  Control  of  this  fundamental 
numerical  instability  is  of  fundamental  importance  in  the  successful  implementation  of 
the  differential  inversion  algorithm. 


3  NUMERICAL  IMPLEMENTATION  OF  THE  DIF¬ 
FERENTIAL  INVERSION  ALGORITHM 

Performing  a  differential  inversion  for  the  determination  of  temperature  as  a  function 
of  altitude,  utilizing  the  infinite  series  derived  in  the  previous  section,  requires  the 
computation  of  the  inversion  coefficients,  A,  ,  and  the  derivatives  (out  to  nth  order)  of 
the  radiance  with  respect  to  f  =  —  ln(p). 

The  derivatives  of  the  radiance  were  estimated  by  a  centered  finite-difference 
technique  (c/.  Chapra  and  Canale  1985).  This  was  carried  out  by  constructing  a  one 
dimensional  grid  around  each  point  for  which  the  derivative  was  to  be  evaluated,  and 
then  applying  the  finite  differences  to  those  points.  To  lowest  order,  the  truncation 
error  of  this  technique  falls  off  as  0(h2),  where  h  is  the  spacing  between  grid  points. 

The  lowest  order  version  of  the  finite  difference  requires  a  total  of  four  grid  points 
in  order  to  obtain  derivatives  to  fourth  order  at  a  single  data  point.  Therefore,  in  the 
particular  problem  addressed  here,  an  atmospheric  radiance  profile  is  decomposed  into 
an  51V  x  2  matrix  consisting  of  points  corresponding  to  ln(p)  and  radiance  over  N 
vertically  distributed  layers  of  the  atmosphere,  and  four  grid  points  symmetrically 
located  about  each  of  the  N  layers  each  also  consisting  of  a  value  of  both  radiance  and 
ln(p). 

The  grid  point  spacing  in  ln(p)  was  initially  taken  as  constant  for  the  synthetic 
data  studies.  However  for  real  data  this  is  not  always  feasible.  Therefore,  another 
variation  in  the  centered  finite-difference  technique  that  accomodates  variable  grid 
point  spacing  must  be  used. 


Pressure  (nb) 


US  STANDARD  ATMOSPHERE:  n=0.6 


Radiance  (erg/s/sq .cw/str ) 

Figure  1:  Synthetic  radiance  profile  computed  for  a  generalized  exponential  weight 
function  with  m  =  0.6.  Radiances  were  calculated  using  the  U.S.  Standard  Atmo¬ 
sphere. 


4  APPLICATION  TO  SYNTHETIC  DATA 


In  order  to  ascertain  the  strengths  and  weaknesses  of  the  King  inversion  equation,  the 
implementation  described  above  was  applied  to  synthetic  radiance  data.  An  evaluation 
of  the  impact  of  truncation  of  Equation  (8)  and  of  numerical  errors  on  the  accuracy 
and  stability  of  the  differential  inversion  technique  is  made  possible  by  examination 
of  the  differential  inversion  independently  of  the  presence  of  experimental  noise  in  the 
data. 

The  radiance  profile  was  computed  from  the  transfer  equation  using  the  gener¬ 
alized  exponential  in  equation  (9)  as  a  weighting  function.  Atmospheric  temperature 
and  pressure  information  specified  as  the  U.S.  Standard  Atmosphere  (Champion  et 
al.  1985).  A  constant  wavenumber  was  used  throughout  the  synthetic  data  study.  In 
order  to  account  for  the  contribution  of  the  earth’s  surface  to  the  upwelling  radiance 
a  surface  term, 

R0  =  D0[l-  W(pfP)dp/p] 

was  included.  Figure  1  shows  the  output  radiance  profile  from  such  a  synthetic  radiance 
profile  calculation. 


In  order  to  accommodate  points  in  the  atmosphere  near  the  surface  of  the  earth, 
for  which  some  of  the  associated  grid  points  may  actually  fall  below  the  surface  of  the 
earth,  the  surface  term  allows  the  radiancp  profile  to  be  extrapolated  by  treating  the 
earth  as  a  semi-infinite,  isothermal  component  of  the  atmosphere. 

The  accuracy  of  the  synthetic  radiance  profile  is  limited  by  numerical  error  in 
the  computations.  All  calculations  were  performed  in  double-precision.  Initially,  the 
transfer  equation  was  integrated  using  Simpson’s  rule;  however,  gross  round-off  errors 
led  to  the  adoption  of  Romberg’s  method  of  integration  (Press  et  al.  1986).  Again,  the 
error  level  was  much  higher  than  expected,  but  nonetheless,  orders  of  magnitude  better 
than  was  obtained  with  Simpson’s  method.  The  error,  primarily  roundoff  from  the 
integration,  is  on  the  order  of  ±0.01  ergs/sec/cm2/ster/cm_1  for  the  synthetic  radiance 
profile.  This  number  is  typical,  but  does  not  apply  to  all  weight  functions.  Later,  in 
the  discussion  it  will  be  shown  that  generalized  exponential  weighting  functions  with 
m  =  0.1  and  m  =  1.0  produce  different  inversion  accuracies  which  result  from  the 
propagation  of  the  error  discussed  here. 

The  radiance  profile  is  fed  directly  into  the  inversion  routine  where  the  derivatives 
are  obtained.  Values  of  m  of  0.1,  0.5,  and  1.0  for  generalized  exponential  weight 
functions  were  used  to  cover  a  range  of  possible  atmospheric  weight  functions.  Inversion 
coefficients  for  these  m  values  were  taken  from  King  (1985). 

The  comparison  of  the  resulting  temperatures  with  those  arising  directly  from 
the  U.S.  Standard  Atmosphere  are  expressed  as  a  simple  difference 

AT  —  T0  —  Tu 

and  an  rms  fluctuation  of  AT  about  zero. 

ATrm,  =  JE(AT)*/N. 

The  spacing  between  points  in  the  finite  difference  grid  h,  proved  to  be  a  parameter  to 
which  the  inversion  is  quite  sensitive.  Figure  2  illustrates  the  accuracy  of  the  inversion 
relative  to  the  stepsize  h.  The  graph  implies  that  the  grid  about  a  single  point  must  be 
large  on  the  scale  of  the  entire  atmosphere  in  order  to  minimize  the  rms  fluctuations 
of  AT  about  zero.  The  problem  a'  es  from  the  fact  that  the  synthetic  radiance  profile 
is  not  error  free.  To  illustrate,  define 

g{xi,h)  =  (/(xl+1)  -  /(x,-j))/2h 

where  f(x)  is  a  differentiable  (order  n)  function,  and  Xi+\  and  x,_i  are  points  equally 
spaced  on  either  side  of  a  point  x,.  For  brevity,  /l  +  1  =  /(x,+i),  /,_ i  =  /(x,_ i),  and 
fi  =  f{Xi). 

Define  a  function  K{x),  which  is  constructed  from  some  point-by-point  represen¬ 
tation  of  /(x)  (».e.  data),  then 
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Figure  2:  Calculation  of  total  numerical  error  on  the  basis  of  synthetic  radiance  calcu¬ 
lations  as  a  function  of  stepsize  h  for  numerical  differentiation.  For  stepsizes  less  than 
approximately  0.35  the  total  numerical  error  begins  to  diverge. 


where  6,  is  the  degradation  contained  in  the  construction  of  K{x),  and  is  defined  as 

6i  =  \K{xi)-f{xi)\. 


Define  a  function 


G(h)  =  (K(xi+1)  -  K[Xi-i))/2h 

where  it+1  and  x,_j  are  as  defined  above.  This  may  be  rewritten  as 


G{h)  =  g(/i)  ±  a, 


with 


<  a  <  a„ 


' min  ^  v  ^  u  max > 


where 


a  =  |<5,+1  ±  $,_,|/2/i, 
°ma  x  =  (|<5.  +  l|  +  |^i-l|)/2/l 


and 


=  0. 


Thus, 


Hrn  G{h)  =  df  jdx  ±o. 


It  can  be  seen  that  a  may  become  very  large  in  the  limit  h  —►  0,  because  the  6's  can 
be  nearly  independent  of  h.  Consider  a  specific  example,  let  all  the  £’s  equal  0.01, 
then  Table  3  gives  values  of  omaz  for  estimates  of  derivatives  of  K(x)  (i.e.  f(: r))  at 
h  =  0.001,  0.01,  0.1,  and  1.0. 


Figures  3  through  5  show  the  inversion  results  on  synthetic  data  with  various  m 
values  for  the  generalized  exponential  weighting  functions.  Point  spacing  for  each  case 
is  h  =  0.3-5  for  m  =  0.1,  h  =  0.70  for  m  —  0.5,  and  h  =  1.1  for  m  =  1.0.  Note  that  the 
fluctuation  about  zero  increases  at  high  pressure.  With  the  stepsizes  indicated  above, 
the  remaining  error  in  the  inversion  stems,  most  likely,  from  the  truncation  of  Equation 
(8),  after  the  fifth  term. 


The  results  just  given  indicate  that  noise  inherent  in  data  will  be  critical  to  the 
application  of  this  inversion  technique.  In  real  data,  one  typically  does  not  have  the 


Table  3:  Total  error  in  the  estimation  of  derivatives  as  a  function  of  stepsize  h. 
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Figure  3:  Errors  in  differential  inversion  of  synthetic  radiance  data  as  a  function  of 
pressure  with  h  =  0.35  and  m  =  0.1. 
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Figure  4:  Errors  in  differential  inversion  of  synthetic  radiance  data  as  a  function  of 
pressure  with  h  =  0.70  and  m  —  0.5. 


luxury  of  positioning  data  points  for  optimization  of  a  single  application  of  that  data. 
In  order  to  gain  some  insight  into  the  behavior  of  King’s  equation  in  the  presence  of 
noise,  random  gaussian  noise  with  selectable  amplitude  was  added  to  the  synthetic 
data. 


Radiance  profiles  with  grid  spacings  of  h  —  0.35  and  h  =  0.9  were  generated 
using  a  generalized  exponential  weighting  function  with  rn  =  0.6.  The  results  of  King's 
differential  inversion  technique  are  given  in  Figures  G  and  7.  Again,  the  critical  role  of 
the  grid  point  spacing,  h,  is  apparent  upon  comparison  of  these  two  curves. 
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Figure  6:  Accuracy  of  differential  inversion  in  the  presence  of  Gaussian  distributed 
random  noise  in  the  radiance  data.  Values  of  m  =  O.G  and  h  =  0.35  were  used. 


ACCURACV  TRADE-OFF 


2.0  4.0  6.0  8.0  10.0 

Random  Noise  (x) 


Figure  7:  Accuracy  of  differential  inversion  in  the  presence  of  Gaussian  distributed 
random  noise  in  the  radiance  data.  Values  of  m  =  0.6  and  h  =  0.95  were  used. 


5  APPLICATION  TO  MEASURED  RADIANCES 


5.1  General  Discussion 


A  basic  purpose  of  this  work  was  to  apply  the  differential  inversion  technique  to  mea¬ 
sured  atmospheric  radiance  data.  This  particular  aspect  of  the  work  had  the  following 
objectives: 


1.  To  define  a  practical  methodology  for  application  of  the  algorithm  to  atmospheric 
measurements. 

2.  To  evaluate  the  performance  of  the  algorithm  in  temperature  inferencing  from 
measurements.  Specifically,  we  sought  the  sensitivity  of  the  algorithm  to  such 
factors  as  spectral  resolution  of  the  measurements,  the  properties  of  atmospheric 
weight  functions,  and  noise  in  the  radiance  measurements. 

3.  To  provide  insight  into  the  kinds  of  atmospheric  radiance  measurements  that 
would  support  operational  temperature  inferencing  by  the  differential  inversion 
technique. 


Two  sources  of  measured  data  were  identified  for  this  application.  The  first  is 
the  Air  Force  Geophysics  Laboratory’s  Stratospheric  Cryogenic  Interferometric  Ex¬ 
periment  (SCRIBE).  The  second  is  the  NOAA  TIROS  Operational  Vertical  Sounder 
(TOVS).  In  section  5.2  we  discuss  the  application  of  the  algorithm  to  the  SCRIBE 
data.  Section  5.3  discusses  the  application  to  the  TOVS  data. 


5.2  The  SCRIBE  Application 


The  SCRIBE  program  was  established  by  the  Air  Force  Geophysics  Laboratory  for  the 
primary  purpose  of  measuring  straospheric  emisssions  at  high  spectral  resolution  from 
balloon  altitudes  (Vanasse  1981;  Gallery  et  al.  1987).  Two  flights  of  the  SCRIBE  (23 
Oct  83  and  5  Jul  84)  at  Holloman  AFB,  NM,  were  identified  as  potential  data  sources 
for  this  study. 

The  SCRIBE  is  a  cryogenically  cooled  Fourier  transform  spectrometer  configured 
to  provide  radiance  data  for  the  spectral  interval  from  500  cm-1  to  2500  cm-1.  The 
spectral  resolution  of  the  measurements  is  about  0.06  cm-1.  The  instrument  may  be 
stepped  through  four  different  viewing  directions  ranging  from  7.5°  above  the  horizon  to 
the  nadir.  Spectral  scans  are  repeated  at  approximately  30  second  intervals  (Murcray 
et  al.  1985,  1984;  Vanasse  1981).  Data  are  transmitted  from  the  balloon  to  a  ground 
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recording  station  during  the  balloon  flight.  Reduction  of  the  digital  interferograms  into 
radiance  spectra  requires  several  processing  steps.  This  processing  was  accomplished 
at  the  University  of  Massachusetts  for  the  data  which  we  used  in  our  application  (Sakai 
1985). 

Gallery  et  al.  (1987)  carefully  examined  the  October  1983  data  for  radiometric 
and  spectroscopic  calibration,  noise  level  and  comparison  of  the  spectral  radiance  dis¬ 
tributions  to  expected  distributions  calculated  from  the  AFGL  FASCOD2  model.  The 
rms  noise  level  of  these  data  has  been  estimated  to  be  of  the  order  of  approximately 
2  x  10~7  watts/(cm  2  sr  cm-1).  In  addition,  the  spectral  registration  of  the  data  is 
known  to  be  displaced  from  its  true  position  by  a  constant  multiplicative  factor.  Based 
on  these  results  we  identified  two  downlooking  scans  from  the  23  October  83  flight  that 
were  suitable  for  application  of  the  differential  inversion  algorithm.  The  two  scans  were 
taken  about  three  minutes  apart.  Table  4  presents  pertinent  details  of  the  two  scans. 
In  order  to  improve  the  spectral  registration  of  the  data,  we  have  analyzed  the  data 
in  reference  to  the  very  accurately  known  positions  of  the  strong  absorption  lines  of 
carbon  dioxide  in  the  spectral  region  from  600  to  1200  cm-1. 


Table  4:  Principal  parameters  describing  the  two  scans  of  SCRIBE  data  from  the  23 
October  1983  flight  that  were  used  to  test  the  application  of  the  differential  inver¬ 
sion  algorithm.  The  two  scans  are  identified  as  SP830V3  and  SP830V1,  respectively. 
(Gallery  et  al.  1987) 


Parameter 

Scan  SP830V3 

Scan  SP830W1 

Minutes  after  launch 

100.8 

104.0 

Altitude  above  mean  sea 
level  (km) 

30.36 

30.67 

Zenith  view  angle 
(degrees) 

180.0 

180.0 

Radiosonde  temperature  at 
flight  altitude  (K) 

232 

233 

667  cm-1  brightness 
temperature  (K) 

225 

225 

800  cm-1  brightness 
temperature  (k) 

270 

272 

Figure  8  shows  a  full  resolution  plot  of  the  SCRIBE  radiances  from  ScanSP830V3 
in  the  spectral  region  at  675  —  681  cm-1.  In  this  opaque  portion  of  the  spectrum, 
the  radiance  profile  represents  the  temperature  structure  immediately  below  balloon 
altitude,  where  the  temperature  increased  with  height  (see  Figure  10).  Here  the  peaks 
in  the  radiance  profile  represent  peak  absorption  positions  corresponding  to  the  strong 
lines.  The  positions  of  the  four  strongest  absorption  lines  in  this  region  are  represented 
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Figure  8:  Radiance  spectrum  from  SCRIBE  Scan  830V3  on  23  July  1983.  The  vertical 
lines  indicate  positions  of  strong  carbon  dioxide  lines  at  676.0196,  677.6010,  679.1856, 
and  680.7733  cm-1  respectively.  Radiance  units  are  erg/crn2sec  ster  cm-1  scaled  by  a 
factor  of  106. 


by  the  vertical  lines. 

Figure  9  shows  a  similar  plot  in  the  spectral  region  at  705  —  712  cm-1.  II  re, 
the  atmosphere  is  more  transparent  and  the  radiant  energy  originates  in  a  layer  where 
the  temperature  decreases  with  height  above  ground.  Hence,  the  strong  line  positions 
are  located  in  the  radiance  minima  in  the  spectrum.  The  strong  line  positions  are  also 
represented  on  this  plot  by  vertical  lines. 


► 


The  extrema  in  the  radiance  spectrum  are  shifted  from  line  positions  by  about 
0.3  cm-1.  We  computed  a  correction  factor  K  for  both  scans  such  that  v'  —  Ku,  where 
v  is  the  reported  spectral  position  of  each  data  point  and  v'  is  the  corrected  spectral 
position.  The  value  of  K  for  Scan  830V3  was  1.00049  and  1.00050  for  the  second  scan, 
indicating  consistency  in  its  value. 

Fluctuations  of  K  about  its  mean  correspond  to  a  variation  in  corrected  spectral 
position  of  about  0.03  to  0.05  cm-1.  This  is  the  same  order  of  magnitude  as  the  resolution 
of  the  data.  Hence,  variations  in  the  individual  estimates  of  fC  may  be  ascribed  largely 
to  resolution  limitations  in  the  measured  radiances. 
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Figure  9:  Radiance  spectrum  from  Scan  830V3  in  the  spectral  region  from  705  to  712 
cm-1.  The  vertical  lines  indicate  positions  of  strong  carbon  dioxide  lines  at  706.5763, 
708.2121,  709.8503,  and  711.4909  cm-1  respectively. 
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Selection  of  a  suitable  set  of  spectral  bands  required  for  such  a  radiance  profile  was 
not  a  straightforward  process.  We  had  no  a  priori  knowledge  of  the  specific  distribution 
of  weight  functions  for  a  specified  set  of  bands.  Furthermore,  we  wanted  to  generate  band 
radiance  sets  at  several  different  bandwidths  in  order  to  evaluate  the  effects  of  spectral 
resolution  on  algorithm  performance. 

5.2.1  FASCODE  Calculations 

The  AFGL  FASCODE  Model  was  used  to  generate  high  resolution  transmittance  profiles 
for  the  atmospheric  column  from  balloon  altitude  to  the  earth’s  surface  for  the  period 
of  the  23  Oct  83  flight  of  the  SCRIBE.  The  FASCODE  model  is  described  in  detail 
in  Clough  et  al.  (1986),  Clough  et  al.  (1981)  and  Smith  et  al.  (1978).  FASCODE  is 
a  highly  developed  program  for  fast,  efficient  computation  of  accurate,  high- resolution 
atmospheric  transmission  profiles.  The  model  assumes  a  layered  atmosphere.  The  calcu¬ 
lations  are  based  on  an  optimum  sampling  at  each  layer  from  the  AFGL  line  compilation 
(Rothman  1981). 

We  input  to  the  model  the  tern  .erature  and  moisture  profiles  reported  in  a  Hol¬ 
loman  radiosonde  observation  taken  at  the  approximate  time  of  the  balloon  flight.  The 
profile  is  shown  in  Figure  10.  Only  the  temperature  and  the  moisture  profiles  were  spec¬ 
ified  in  the  input  profile.  US  Standard  Atmosphere  default  profiles  were  used  for  the 
other  constituent  gases  represented  in  the  model. 

In  order  to  achieve  as  much  detail  as  possible  in  the  transmittance  profile  we  gen¬ 
erated  profiles  at  60  levels,  the  maximum  permitted  by  FASCODE,  between  balloon 
altitude  and  the  ground.  This  resulted  in  0.5  km  vertical  spacing  between  transmittance 
levels. 

The  transmittance  profiles  used  in  the  radiance  band  selection  were  constructed 
from  the  0.006  cm-1  transmittances  by  using  the  FASCODE  filter  option.  We  con¬ 
structed  square  filter  functions  with  101  data  points  per  filter.  Five  sets  of  filters  were 
used  to  generate  the  five  transmittance  files  at  the  five  different  values  of  spectral  res¬ 
olution.  Each  transmission  profile  was  differentiated  to  obtain  the  weight  functions 
(dr/ din p)  and  to  locate  the  pressure  level  of  the  weight  function  peak. 

5.2.2  Calculation  of  Weight  Functions  and  Inversion  Coefficients 

Weight  functions  were  calculated  from  the  transmittances  calculated  from  FASCODE  by 
numerical  differentiation.  A  centered  divided  difference  calculation  modified  for  unevenly 
spaced  data  was  used  to  obtain  the  numerical  derivatives.  It  was  determined  that  inac¬ 
curacies  in  the  FASCODE  transmittances  prevented  calculation  of  physically  meaningful 
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Figure  10:  The  atmospheric  temperature  profile  from  the  Holloman  rawinsonde  obser¬ 
vation  of  23  July  1983,  The  approximate  time  of  the  rawindonde  was  0700  MST  (1400 
GMT). 
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weight  functions  throughout  the  range  of  atmospheric  pressures  observed  with  SCRIBE. 
The  problem  is  especially  severe  for  the  small  transmittances  at  the  low  pressure  end  of 
the  FASCODE  pressure  grid  (see  Appendix  A). 

These  inaccuracies  in  FASCODE  transmittances  cause  oscillations  in  the  calculated 
weight  functions  with  a  period  of  two  pressure  grid  intervals.  The  amplitude  of  the 
oscillations  may  be  as  great  as  the  value  of  the  nominal  peak  of  the  weight  function. 
Since  we  expect  that  the  weight  functions  must  be  smooth  functions  of  the  pressure, 
the  weight  functions  obtained  by  numerical  differentiation  of  FASCODE  transmittances 
are  unsuitable  for  use  in  differential  inversion.  The  problem  is  especially  severe  for 
the  evaluation  of  the  inversion  coefficients,  A,-,  which  involve  determination  of  the  ith 
moment  of  the  weight  function  with  respect  to  pressure.  Numerical  errors  at  the  edge  of 
the  pressure  grid  particularly  corrupt  the  evaluation  of  A^  for  large  values  of  i. 

Since  these  oscillations  of  the  weight  functions  calculated  from  FASCODE  trans¬ 
mittances  are  nonphysical,  we  are  justified  in  using  a  fitting  procedure  based  upon  the 
physically  reasonable  portion  of  the  weight  function  curves  obtained  by  numerical  dif¬ 
ferentiation.  We  fitted  generalized  exponential  weight  functions  (see  eq.  (9))  to  these 
FASCODE  derived  weight  functions.  Typically  five  points  from  the  low  pressure  end  of 
the  FASCODE  pressure  grid  were  discarded  in  the  process  of  constructing  the  general¬ 
ized  exponential  fit.  Given  the  m  value  resulting  from  the  generalized  exponential  fit, 
values  of  Aj  may  be  obtained  either  numerically  (using  a  Romberg  integration  algorithm) 
or  on  the  basis  of  analytic  formulae. 

5.2.3  Description  of  SCRIBE  Band  Radiance  Sets 

Two  spectral  scans  from  the  23  October  83  flight  were  selected  for  application  of  the 
differential  inversion  algorithm.  Their  identifications  are  SP830V3  and  SP830W1  re¬ 
spectively.  Table  5  summarizes  the  band  radiance  sets  that  were  constructed  from  these 
two  scans.  The  radiance  sets  were  constructed  by  convolution  of  rectangular  filter  func¬ 
tions  with  the  full  resolution  spectra.  Each  set  contained  anywhere  from  7  to  70  radiances 
depending  on  the  spectral  resolution  aM  the  number  of  radiances  assigned  for  evaluation 
of  derivatives  at  the  p  level.  The  number  of  p  levels  in  each  set  varied  from  1  to  10. 

The  sets  are  ordered  as  follows: 

1.  by  the  two  scans,  SP830V3  and  SP830W1. 

2.  by  spectral  resolution;  the  values  of  resolution  are  indicated  for  each  set. 

3.  by  two  differing  values  of  nominal  vertical  spacing,  i.e.,  by  the  coarseness  of  the 

finite  difference  gird  used  in  numerical  evaluation  of  the  derivatives. 


Table  5:  Summary  of  SCRIBE  Band  Radiance  Sets 


File  name 


BNRD21 

BNRD22 

BNRD23 

BNRD24 

BNRD25 


CNRD21 

CNRD22 

CNRD23 

CNRD24 

CNRD25 


DNRD21 

DNRD22 

DNRD23 

DNRD24 

DNRD25 

ENRD21 

ENRD22 

ENRD23 

ENRD24 

ENRD25 


FNRD21 

FNRD22 

FNRD23 

FNRD24 

FNRD25 


GNRD21 

GNRD22 

GNRD23 

GNRD24 

GNRD25 


Spectral 

Resolution 


Nominal  h 


0.20  to  0.25 


Number  of 
radiances 
per  p 


Number 
p  levels 


0.35  to  0.40 


0.35  to  0.40 


0.20  to  0.25 


0.35  to  0.40 


0.35  to  0.40 


Table  6:  Comparison  of  corresponding  radiances  between  two  scans  of  the  SCRIBE. 
Columns  4  and  5  are  the  radiances  for  the  identified  scans  in  units  of  erg/m  sec  sr  cm-1. 


p 

(mb) 

Bandwidth 

cm-1 

Band  Center 
cm-1 

SP830V3 

SP830W1 

35.510 

0.1 

682.55 

37.569 

36.662 

41.538 

681.05 

36.089 

35.480 

57.263 

672.55 

49.835 

54.627 

79.283 

689.05 

36.065 

37.669 

101.480 

684.55 

38.112 

37.047 

138.540 

701.55 

44.190 

44.281 

192.060 

701.05 

46.003 

43.264 

261.400 

702.55 

54.850 

60.027 

302.820 

712.05 

66.209 

71.770 

490.160 

713.55 

59.914 

60.202 

152.190 

2.0 

692.00 

38.685 

40.212 

192.740 

697.50 

42.042 

41.725 

261.400 

701.00 

45.837 

45.195 

325.300 

704.00 

50.720 

51.250 

418.810 

707.50 

56.739 

58.390 

195.380 

10.0 

695.00 

42.125 

42.919 

230.010 

701.50 

47.200 

47.524 

265.040 

701.00 

46.779 

47.253 

302.820 

702.50 

48.618 

49.217 

4.  by  the  number  of  radiance  levels  assigned  to  each  p  for  numerical  differentiation. 
Seven  radiances  were  used  in  general.  By  reducing  to  five,  we  were  able  to  increase 
by  two  the  number  of  p  levels  in  each  case. 

Table  6  shows  a  comparison  of  the  two  corresponding  radiance  sets  taken  from  the 
two  successive  scans,  for  several  specific  bands  in  the  overall  set.  Note  that  each  pair  in 
this  table  represents  two  successive  measurements  about  three  minutes  apart  in  a  specific 
waveband.  Differences  between  radiances  within  each  pair  are  an  indication  of  radiance 
noise,  since  the  scans  should  see  nearly  the  same  atmospheric  thermal  profile.  We  see 
that  at  0.1  cm-1  resolution,  the  nominal  variation  is  of  the  same  order  as  the  noise 
estimate  reported  by  Gallery  et  al.  (1987),  i.e.,  of  order  2  (erg/cm2  sr  sec  cm-1).  As  the 
bandwidth  increases,  we  note  that  the  differences  within  radiance  pairs  are  smaller  as  a 
result  of  the  averaging  effect  of  the  wider  bands. 


We  have  seen  in  the  algorithm  definition  (Section  3)  and  the  synthetic  radiance 
inversions  (Section  4)  that  a  well  behaved  radiance  profile  of  radiance  vs.  p  is  required 
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Figure  11:  Radiance  values  us.  p  at  four  values  of  spectral  resolution  for  the  SCRIBE 
radiance  data  selected  for  the  inversions. 


for  evaluation  of  the  derivatives.  In  order  to  evaluate  the  utility  of  the  radiances  we 
plotted  the  radiances  within  each  of  the  first  five  files  in  Table  5.  Figure  11  shows  the 
plots.  We  see  that  a  considerable  scatter  occurs  for  the  fine  resolution  data;  at  the  same 
time  we  see  a  well-behaved  pattern  in  the  radiances  with  a  wider  bandwidth. 


We  ascribe  the  scatter  in  the  radiance  data  plots  to  the  causes  listed  below.  The 
relative  importance  of  each  factor  cannot  be  evaluated  quantitatively  at  presend.  How¬ 
ever,  we  have  subjectively  ordered  the  factors  below  (most  important  first)  based  on  our 
understanding  of  the  data  and  our  experience  with  the  FASCODE  model.  These  are  our 
best  estimate  of  the  related  factors. 


1.  Noise  in  the  radiance  measurements.  This  is  undoubtedly  the  primary  factor.  The 
other  indicators  of  noise  level  in  the  SCRIBE  data  indicate  that  a  large  part  of  the 
scatter  can  be  ascribed  to  the  inherent  noise. 
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Table  7:  Results  of  temperature  retrievals  by  application  of  the  differential  inversion 
algorithm  to  the  SCRIBE  radiances  with  bandwidth  equal  to  5.0  cm-1.  The  symbol 
rid  denotes  the  number  of  derivatives  used  in  the  retrieval.  Temperatures  are  expressed 
in  degrees  Kelvin.  Pressure  levels  (p)  are  in  millibars. 


Scan 

rid 

P 

Computed  T 

T(Raob) 

AT 

SP830V1 

4 

192 

234 

218 

+  16 

242 

209 

226 

-17 

3 

192 

236 

218 

+  18 

242 

213 

226 

-13 

SP830W1 

4 

192 

235 

218 

+  17 

242 

202 

226 

-14 

3 

192 

233 

218 

+  15 

242 

188 

226 

-26 

2.  Uncertainty  in  spectral  registration  of  the  radiance  data.  It  is  apparent  from  the 
fine  structure  of  the  radiance  data  in  Figures  11  that  very  slight  errors  in  the 
spectral  position  of  the  radiances  would  produce  substantial  scatter  in  the  narrow 
band  radiances. 

3.  Uncertainty  in  determination  of  p  (or  p),  the  pressure  level  of  the  peak  in  weight 
functions.  The  60  level  limit  in  FASCOD  as  well  as  numerical  techniques  used  to 
locate  p  limit  the  accuracy  of  this  parameter. 


5.2.4  SCRIBE  Radiance  Inversion  Results 

The  results  of  the  differential  inversion  of  the  SCRIBE  data  are  given  in  Table  7.  Owing 
to  the  small  number  of  independent  p  values,  only  two  points  of  the  temperature  profile 
were  obtained.  We  performed  inversions  for  each  of  the  two  SCRIBE  scans  with  the 
differential  inversion  series  truncated  at  third  and  fourth  order,  thus  demonstrating  that 
the  SCRIBE  inversions  are  insensitive  to  the  order  at  which  the  series  is  truncated.  It 
can  be  seen  that  errors  in  the  determination  of  temperatures  were  ~  10  —  20  K,  compared 
to  radiosonde  observations.  The  inversion  of  the  SCRIBE  data  were  carried  out  with 
an  effective  stepsize,  h  «  0.3,  and  so  the  errors  obtained  are  not  inconsistent  with  the 
results  of  the  synthetic  inversions  (see  Table  3).  We  note  that  with  the  small  number  of 
independent  data  points  (4)  given  here  and  their  large  errors  compared  with  radiosonde 
observations,  that  there  is  no  significance  to  the  trend  of  temperature  with  height  shown 
by  differential  inversion  of  this  data. 

The  inversions  of  SCRIBE  data  shown  here  are  not  conclusive  for  demonstrating  the 
efficacy  of  the  differential  inversion  algorithm.  The  level  of  accuracy  of  these  results  illus- 
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trates  problems  for  the  differential  inversion  algorithm  arising  from  the  calibration  and 
noise  level  of  the  SCRIBE  data,  fundamental  difficulties  with  finite-difference  schemes 
in  the  presence  of  noise,  and  the  approximation  of  SCRIBE  weight  functions  by  gen¬ 
eralized  exponentials.  Improved  techniques  for  differential  inversions,  currently  under 
development,  will  control  a  large  portion  of  these  errors. 

5.3  The  TOYS  Application 

5.3.1  Description  of  the  TOVS 

The  TIROS  Operational  Vertical  Sounder  (TOVS)  consists  of  three  instruments:  (1)  the 
High  Resolution  Infrared  Sounder  (HIRS-2),  (2)  the  Stratospheric  Sounding  Unit  (SSU), 
and  (3)  the  Microwave  Sounding  Unit  (MSU).  In  total  the  TOVS  measures  radiances 
in  27  channels.  The  band  limits  of  the  various  channels  and  information  regarding  the 
bandwidths  and  peak  levels  of  the  weight  functions  are  given  in  the  NOAA  Polar  Orbiter 
Data  User’s  Guide  (NOAA  1983)  and  in  Smith  et  al.  (1979).  Of  the  27  TOVS  channels, 
12  channels  of  the  HIRS-2  (1-7  and  13-17),  the  three  SSU  channels,  and  three  of  the  four 
MSU  channels  respond  primarily  to  changes  in  the  atmospheric  temperature  profile. 

The  infrared  channels  of  the  TOVS  cover  fairly  broad  spectral  bands.  For  Channel 
1,  which  has  the  smallest  bandwidth  among  the  infrared  channels,  the  half-power  band¬ 
width  is  about  4  cm-1  and  the  one  tenth-power  bandwidth  is  about  9  cm-1.  The  other 
infrared  channels  have  half-power  bandwidths  of  approximately  20  cm-1  and  one  tenth- 
power  bandwidths  of  approximately  30  cm-1.  Note  that  in  comparison  to  the  SCRIBE 
data  discussed  in  Section  5.2,  the  TOVS  infrared  data  have  comparatively  coarse  spectral 
resolution. 


5.3.2  Selection  of  Data 

TOVS  data  were  provided  by  the  National  Environmental  Satellite  Data  and  Informa¬ 
tion  Service  (NESDIS)  from  the  operational  data  processed  by  NESDIS  in  routine  daily 
operations.  The  particular  radiance  data  which  we  used  in  this  study  were  taken  from 
a  data  sample  collected  during  the  period  16-18  October  1987.  The  sample  consisted 
of  15  clear  column  radiance  sets  from  each  of  three  latitude  bands.  The  radiance  data 
from  each  latitude  band  were  accompanied  by  a  set  of  transmittance  profiles  that  wure 
computed  by  NESDIS  from  the  average  temperature-moisture  profile  in  that  latitude 
band.  These  transmittances  were  computed  by  the  method  developed  by  Weinreb  et  al. 
(1981). 
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5.3.3  Differential  Inversion  of  TOYS  data 


The  transmittances  provided  by  NESDIS  on  a  grid  of  40  atmospheric  levels  were  numeri¬ 
cally  differentiated  to  obtain  weight  functions.  Computation  of  inversion  coefficients,  A, , 
for  large  values  of  t,  is  still  ill-conditioned  since  the  evaluation  of  high  moments  of  the 
weight  function  on  a  coarse  grid  is  required.  Therefore,  we  fit  generalized  exponential 
weight  functions  to  the  weight  functions  found  from  NESDIS  transmittances  and  used 
those  fitted  weight  functions  to  calculate  A {  values.  The  generalized  exponential  fits 
were  obtained  by  first  determining  the  maximum  value  of  the  empirical  weight  function 
based  on  numerical  differentiation  (centered  divided  difference)  of  the  NESDIS  trans¬ 
mittances  and  fixing  this  value  as  p  in  V  =  p/p  as  in  eq.  (9).  Then  an  m  value  was 
determined  which  minimized  the  rms  error  between  the  generalized  exponential  weight 
function  on  the  grid  of  the  NESDIS  pressure  levels.  The  fitting  procedure  exploits  the 
character  of  the  generalized  exponential  function  as  a  one-parameter  family  of  functions 
(i.e.  determined  by  m  values). 

The  differential  inversion  of  TOVS  data  produced  results  significantly  better  than 
those  obtained  for  SCRIBE.  The  procedure  for  the  differential  inversion  of  TOVS  data 
was  identical  to  that  for  the  SCRIBE.  Pressure  levels  were  selected  to  give  the  most 
stable  results  (see  discussion  in  sections  4  and  5.2.4),  with  a  stepsize  h  ft;  0.3  —  0.4.  Table 
8  gives  the  temperature  values  obtained  for  the  2190  —  2360  cm"1,  mid-latitude  TOVS 
data.  Figure  12  follows  by  plotting  inferred  temperatures  against  Raob  temperatures  for 
all  TOVS  data.  The  discrepancy  between  inferred  temperatures  and  Raob  data  is  on  the 
order  of  ±2  —  3  degrees  for  points  in  the  p  grid  selected  in  the  2190  to  2360  cm-1  range 
and  ±15  —  20  degrees  for  points  in  the  668  to  748  cm-1  region. 

For  sake  of  completeness,  we  shall  show  explicitly  the  computations  for  one  of  the 
TOVS  differential  inversions.  The  radiance  \alues  for  the  first  scan  in  Table  8  are  given 
in  Table  9.  For  a  generalized  exponential  weight  function  with  m  =  0.49,  we  obtain  by 
Romberg  integration  differential  inversion  coefficients  of  Ao  =  1.0,  Ai  =  —0.621798,  A2  = 
—0.451558,  and  A3  =  0.00438163.  Evaluating  the  radiance  derivatives  and  multiplying 
by  the  differential  inversion  coefficients,  we  find, 

B(4Q0  mb)  =  0.165079  +  0.120947  -  0.1  3742  -  0.000557984  =  0.11726  erg/s  cm2  str  cm'1  (1) 

using  terms  out  to  third  order.  This  result  suggests  that  convergence  is  rapid,  in  fact,  the 
third  order  term  does  not  contribute  significantly  to  the  derived  temperature.  This  rapid 
convergence  implies  that  application  of  differential  inversion  will  be  computationally 
effective  in  practical  applications. 

We  note  that,  within  the  context  of  very  limited  data,  the  TOVS  results  in  the  668 
to  748  cm-1  region  are  similar  to  the  temperature  differences  obtained  for  the  SCRIBE 
data.  The  difference  in  behavior  between  the  wavelength  regions  suggests  that  at  small 
inverse  wavelength,  small  changes  in  radiance  correspond  to  large  changes  in  tempera- 


Table  8:  Differential  inversion  results  for  TOVS.  The  value  of  p  is  400  millibars. 


Latitude  Longitude  Time 


-82.6 

123.8 
122.1 

120.8 
-171.4 

32.6 

139.4 
-82.6 
154.2 
-79.7 
-106.9 

153.4 
129.1 
104.9 


13:30:22 

11:12:19 

11:13:51 

23:41:12 

18:46:07 

17:42:47 

22:09:22 

00:44:30 

09:09:21 

13:08:17 

14:50:41 

21:35:20 

23:15:31 

12:29:12 


Inferred 

Temperature 

253.27 

255.53 

254.83 

256.00 

255.04 

247.87 

253.57 

254.40 

257.02 


Raob 

Temperature 


252.75 

256.25 


253.61 
253.63 
256.88 

253.62 
256.79 


255.25 
257.16 
257.20 
245.19 

254.25 
252.05 
258.86 
250.55 
255.69 
258.45 
252.95 

258.25 


+  1.52 
-0.72 
-0.42 
-1.16 
-2.16 
+2.68 
-0.68 
+2.35 
-  1.84 
+3.06 
-1.06 
-1.57 
+0.67 
-1.46 


TOUS  INUERSION  RESULTS 


Inferred  Tenperature  (H) 
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Figure  12:  Temperature  values  obtained  from  application  of  the  differential  inversion 
algorithm  to  TOVS  data  are  plotted  against  corresponding  Raob  values.  Note  that  the 
low  wavenumber  temperatures  are  systematically  offset  from  the  Raob  measurements. 
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Table  9:  Radiance  values  for  a  representative  TOVS  scan,  used  for  an  example  differential 
inversion  computation  carried  out  in  the  text.  Radiances  (erg/cm2  sec  sr  cm-1)  have 
been  normalized  to  channel  15. 


Channel 

Wavenumber 

(cm-1) 

P 

(millibars) 

Brightness 
Temp.  (K) 

Radiance  (2238.45 

H 

2190.10 

990.45 

280.360 

1.37306 

2195.10 

1068.75 

269.410 

0.860927 

15 

2238.45 

400.00 

250.230 

0.344447 

16 

2264.95 

175.08 

230.950 

0.117654 

17 

2361.70 

20.00 

239.470 

0.193209 

ture,  thus  temperature  inference  accuracy  will  be  much  more  strongly  affected  by  noise 
than  the  large  inverse  wavelength  region  of  the  radiance  profile. 
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6  Conclusions 


Realistic  numerical  experiments  on  synthetic  radiance  data  have  demonstrated  that  dif¬ 
ferential  inversion  techniques  may  be  successfully  employed  to  obtain  atmospheric  tem¬ 
perature  profiles  from  radiance  data.  Applications  of  the  differential  inversion  technique 
to  SCRIBE  and  TOVS  radiance  data  have  shown  practical  problems  that  can  be  re¬ 
solved  by  improvements  in  the  application  of  the  differential  inversion  technique  and 
in  the  selection  of  instrumentation.  At  least  some  of  these  difficulties  were  overcome 
in  our  inversions  of  TOVS  radiance  data,  and  the  direction  of  development  for  further 
improvements  is  clear. 

When  implemented  in  such  a  fashion  as  to  require  numerical  differentiation  of 
radiance  data,  differential  inversion  is  sensitive  to  the  presence  of  noise  and  to  calibration 
errors  in  the  radiance  data.  These  problems  limit  the  stepsize,  h,  of  the  calculations  in 
logarithmic  pressure  space  so  as  to  control  the  magnitude  of  the  relative  numerical  error. 
In  application  to  SCRIBE  data,  differential  inversion  produced  typical  errors  ~  20  K, 
and  for  TOVS  data  typical  errors  ~  2  K. 

The  problems  outlined  here  are  fundamental  to  any  form  of  the  differential  inversion 
algorithm  dependent  upon  finite-difference  evaluation  of  derivatives.  Improvements  in 
the  application  of  the  differential  inversion  technique  based  upon  insights  gained  from 
radiative  transfer  theory  may  allow  us  to  circumvent  this  difficulty. 


7  Recommendations 


As  a  result  of  our  research  it  is  possible  to  make  several  recommendations  which  will 
allow  more  effective  exploitation  of  the  advantages  of  the  differential  inversion  algorithm, 
particularly  its  computational  efficiency,  in  future  generations  of  atmospheric  sounders. 

1.  A  larger  number  of  channels  is  required  to  evaluate  the  temperature  profile  at  more 
heights  in  the  atmosphere,  particularly  if  finite  difference  schemes  are  employed 
(but  see  item  4  below). 

2.  Accurate  calibration  and  high  spectral  resolution  of  the  atmospheric  sounding  hard¬ 
ware  will  permit  more  accurate  evaluation  of  numerical  derivatives,  with  a  substan¬ 
tial  positive  impact  on  the  accuracy  of  the  differential  inversion  results. 

3.  If  a  new  instrument  is  to  be  designed,  an  investigation  into  the  availability  of 
channels  that  can  be  selected  such  that  narrow  weight  functions  corresponding  to 
desired  p  values  can  be  obtained. 

4.  Algorithmic  techniques  eliminating  the  requirement  for  evaluation  of  numerical 
derivatives  in  the  differential  inversion  technique  should  be  vigorously  pursued. 
These  theoretical  developments  could  have  the  largest  positive  impact  on  the  utility 
of  the  differential  inversion  algorithm. 
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A  SCRIBE  Weight  Functions  from  TOVS  Trans- 
mittances 


Initially  weight  functions  used  for  differential  inversion  were  computed  directly  by  nu¬ 
merical  differentiation  of  a  transmittance  profile  computed  using  FASC0D2.  A  pressure 
grid  of  60  levels  was  used,  which  is  the  finest  grid  supported  by  FASCOD2.  The  weight 
functions  were  computed  by  a  standard  centered  divided  difference  scheme,  with  for¬ 
ward  and  backward  divided  difference  calculations  performed  at  the  edge  of  the  grid,  as 
required. 

In  the  figures,  a  randomly  selected  sample  of  weight  functions  computed  from  FASCOD2 
transmittances  is  shown  both  with  pressure  and  logarithm  (base  10)  of  pressure  as  inde¬ 
pendent  variable.  The  plot  of  weight  function  against  logarithm  of  pressure  is  particularly 
useful  for  illustrating  the  problems  with  the  weight  functions  we  have  found  here,  which 
are  most  important  at  low  pressures.  Although  numerical  inaccuracies  of  the  FASCODE 
transmittances  are  not  apparent  in  plots  of  the  transmittance  function,  tney  are  sig¬ 
nificant  for  the  weight  functions,  because  numerical  differentiation  effectively  amplifies 
relative  numerical  error.  The  result  is  that  the  weight  functions  are  no  longer  smooth, 
but  exhibit  substantial  fluctuations  between  adjacent  points  over  the  low  pressure  range 
of  the  pressure  grid. 

Application  of  the  differential  inversion  algorithm  assumes  that  the  weight  functions 
must  be  reasonably  smooth  functions,  and  in  fact  this  is  an  important  requirement  for 
the  numerical  stability  of  the  inversion  algorithm.  A  portion  of  the  difficulty  may  be 
seen  by  considering  the  evaluation  of  the  A  coefficients  for  the  series  in  equation  (8). 
Determination  of  these  coefficients  requires  the  evaluation  of  moments  of  the  weight 
function.  Applying  equation  (6), 


which  implies 
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and  so  evaluation  of  A*.  involves  computation  of  the  k  —  1th  and  lower  moments  of  the 
weight  function.  For  large  values  of  k,  calculation  of  moments,  particularly  emphasizes 
the  importance  of  any  inaccuracies  in  the  “wings”  of  the  weight  function,  t.e.  at  the  low 
pressure  and  high  pressure  and  high  pressure  ends  of  the  pressure  grid.  Therefore,  the 
calculation  of  A  coefficients  from  FASCODE  transmittances  is  unreliable. 


To  circumvent  this  problem  for  application  of  differential  inversion  to  the  SCRIBE  data, 
we  fitted  generalized  exponential  weight  functions  to  the  weight  functions  obtained  by 
numerical  differentiation  of  the  FASCOD2  transmittances.  Values  of  m  were  determined 
to  three  significant  figures  by  minimization  of  mis  errors  on  the  FASCODE  pressure  grid. 

36 


£ 


Using  this  m  value,  A  coefficients  were  found  for  the  inversion  by  taking  appropriate 
moments  of  the  generalized  exponential  weight  functions. 

Inspection  of  the  graphs  of  the  weight  functions  shows  that  the  problems  with  the  FAS- 
COD2  transmittances  principally  arise  at  the  edge  of  the  pressure  grid  and  might  be 
attributed  to  some  sort  of  edge  effect  in  the  FASCOD2  computations  of  the  transmit¬ 
tances.  The  fluctuations  in  the  weight  functions  are  not  attributable  to  use  of  forward 
or  backward  divided  differences  at  the  edge  of  the  grid  for  computation  of  the  weight 
function  (as  opposed  to  centered  divided  differences)  since  the  fluctuations  occur  up  to 
5  or  6  grid  points  interior  to  the  edge  of  the  pressure  grid.  The  relative  accuracy  of 
computation  of  FASCODE  transmittances  will  have  to  be  increased  by  approximately 
an  order  of  magnitude  in  order  to  achieve  the  degree  of  smoothness  required  to  directly 
determine  A  coefficients  from  FASCODE  transmittances. 

Some  advantage  may  be  gained  by  smoothing  the  FASCODE  transmittances  to  the  actual 
accuracy  of  the  transmittances.  This  smoothing  would  alleviate  some  of  the  problems 
arising  from  the  fluctuations  in  the  transmittance  profiles  at  the  low  pressure  end  of  the 
pressure  grid.  This  same  smoothing  effect  is  probably  inherent  to  the  NOAA  statistical 
transmittance  model. 
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Figure  Al:  Randomly  selected  weight  functions  determined  by  numerical  differentiation 
of  FASCODE  transmittances.  In  (a)  the  independent  variable  is  the  atmospheric  pressure 
in  bars,  whiie  in  (b)  it  is  the  log10  of  the  pressure  in  millibars. 


